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PREFACE 


The  main  objective  of  this  study  is  to  develop  a  two-dimensional 
numerical  model  capable  of  studying  water  and  sediment  movement  and 
geomorphic  changes  in  alluvial  channel  reaches  with  complex  geometries.  A 
description  of  the  governing  equations  for  water  and  sediment  motion  in  two 
dimensions  is  presented  in  Part  2  of  this  report.  Part  3  presents  a 
detailed  description  of  the  numerical  methods  used  to  discretize  the 
governing  equations  so  that  the  solution  can  be  carried  out  with  the  aid  of 
a  digital  computer.  Part  A  of  this  report  presents  results  of  applications 
of  the  model  to  the  prediction  of  a  trench  scour  and  fill,  and  the 
simulation  of  bed  scour  around  a  spur  dike. 
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INTRODUCTION 


The  objective  of  this  study  is  to  develop  a  new  methodology  for 
modeling  the  phenomenon  of  sediment  movement  in  irregular  alluvial 
channels,  scouring  around  obstructions,  etc.  The  basic  physical  principles 
of  conservation  of  mass  and  momentum  are  used  to  describe  the  fluid  flow. 
The  conservation  of  mass  and  semi-empirical  equations  governing  sediment 
particle  movement  are  adopted  to  establish  the  interaction  between  the 
sediment  movement  and  the  fluid  flow.  The  resulting  mathematical  model  is, 
unfortunately,  highly  nonlinear  and  complex.  It  is  impractical,  if  not 
impossible,  to  solve  it  analytically.  Therefore,  the  numerical  methods  of 
finite  element  and  finite  difference  are  used  to  obtain  approximate 
solutions  of  this  model. 

The  application  of  the  finite  element  method  (FEM)  to  model  fluid 
flows  has  progressed  rapidly  in  recent  years  from  the  simplest  linear 
inviscid  fluid  flow  problems  (Martin,  1968;  Argyris  et  al.,  1969)  to  slow 
viscous  flows  (Tong,  1969;  Atkinson  et  al.,  1969;  Oden  and  Sornogyi,  1969), 
and  finally  to  the  solution  of  the  full  Navier-Stokes  equations  (Oden, 
1970;  Skiba,  1970;  Olson,  1972;  Oden  and  Wellford,  1972).  However,  this 
latter  area  represents  an  extremely  large  and  complex  field.  As  such,  the 
research,  although  very  active,  can  only  be  referred  to  as  being  in  its 
beginning  stage  (Olson,  1975).  A  summary  of  its  recent  applications  to 
flows  through  porous  media,  shallow  water  circulation,  and  two-dimensional 
viscous  flows  had  been  presented  by  Connor  and  Brebbia  (1976) .  ^Jiorrie  and 
de  Vries  (1978)  also  presented  an  excellent  survey  of  the  FEM  applications 
in  all  branches  of  fluid  mechanics  with  218  papers  cited.  Readers, 
desiring  to  find  detailed  information  on  the  development  of  FE  Modeling  of 
Fluid  Flows  in  general,  are  referenced  to  these  and  other  similar  papers. 
A  complete  review  on  FE  Modeling  of  Open  Channel  Flows  and  directly  related 
works  is  presented  below. 

A  variational  principle  for  an  ideal  fluid  flow  with  a  free  surface 
under  gravity  was  developed  by  Luke  (1967)  using  potential  function 
formulation.  It  was  modified  using  the  stream  function  formulation  and  the 
different  expressions  of  free  surface  boundary  conditions  by  O'Carroll 
(1975),  and  by  O'Carrol  et  al.  (1976).  O'Carrol  (1975a,  1978)  also  applied 
this  FEM  to  compute  the  flows  by  a  vertically  two-dimensional  model  and 


over  a  spillway,  etc.  Although  he  only  solved  the  problems  without  the 
effects  of  side  walls;  he  contributed  a  great  deal  to  techniques  for 
handling  the  moving  free  surface  at  least  for  the  ideal  fluid  flow.  Using 
the  Galerkin's  approach  of  the  FEM,  Keuning  (1976)  solved  a  straight 
horizontal  channel  with  a  uniform  trapezoidal  cross  section.  Although  the 
problem  is  only  one-dimensional,  the  equations  are  kept  nonlinear  and 
unsteady . 

A  two-dimensional  quasi-linear  FE  Model  for  Open  Channel  Flow  near 
Critical  Conditions  was  reported  recently  by  Katopodes  (1980).  It 
successfully  demonstrated  the  capability  of  FEM  to  simulate  a  supercritical 
floodwave.  The  truly  three-dimensional  finite  element  modeling  of  viscous 
flows  in  an  open  channel  with  and  without  the  existance  of  obstructions  was 
carried  out  by  Alonso  and  Wang  (1978).  Three-dimensional  linear 
hexahedral,  isoparametric  elements  were  used  to  obtain  very  slow  viscous 
laminar  flows  in  open  channels  of  varying  cross-section  and  around  an 
isolated  obstruction.  Although  results  obtained  were  physically  sound  and 
mathematically  reasonable,  the  requirement  of  computer  storage  and 
computing  time  were  prohibitive.  One  of  the  most  effective  alternatives  is 
the  depth-averaging  scheme.  It  has  been  used  primarily  in  the  simulation 
of  currents  and  water  waves  in  lakes,  estuaries  and  shallow  coastal  water. 
Some  typical  contributions  may  be  found  in  Leendertse  (1967),  Nakayama  and 
Romero  (1971),  and  Niemeyer  (1977).  More  recently  the  utilization  of  the 
depth-averaging  models  in  the  finite  element  simulation  of  flows  in  open 
channels  and  rivers  were  reported  by  Thienpont  and  Berlamont  (1980)  and 
Wang,  Su,  and  Alonso  (1980).  Because  the  distribution  of  hydrodynamic 
properties  in  the  vertical  direction  of  a  shallow  water  flow  are  usually 
better  understood,  appropriate  functions  can,  thus,  be  chosen  to  yield 
adequate  approximation  in  this  direction.  Therefore,  the  governing 
differential  equations  can  be  integrated  vertically  from  the  channel  bed  to 
the  free  surface  resulting  in  differential  equations,  containing  vertically 
averaged  properties,  which  are  only  two-dimensional  in  a  horizontal 
reference  plane.  Even  if  the  time  derivatives  are  retained  in  the 
equations  to  model  unsteady  flows,  the  requirement  of  computer  storage  as 
well  as  computing  time  to  simulate  an  open  channel  flow  is  greatly  reduced. 
Besides,  this  approach  not  only  gives  reasonable  results  with  adequate 
accuracy,  but  allows  better  resolution  in  horizontal  directions  by  using 
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the  computer  storage  saved  from  reducing  three-dimensional  to  two- 
dimensional  formulation  to  add  more  nodes  on  a  horizontal  plane. 
Furthermore,  the  computer  code  developed  based  on  this  approach  has  the 
potential  of  wider  acceptance  by  users  with  limited  computing  resources. 
The  simulation  of  sediment  transport  is  discussed  below. 

In  recent  years  one-dimensional,  mathematical  models  of  sediment 
routing,  morphological  transients,  and  sediment  deposition,  etc.  were 
developed  (Cunge  and  Perdreau,  1973;  Mahmood,  1979;  Simons  et  al.,  1975; 
Lopez,  1978).  Although  they  do  not  provide  the  time-varying  configuration 
of  the  sand  bed  in  a  horizontal  plane,  these  models  contribute  a  great  deal 
in  understanding  the  basic  characteristics  of  morphological  transients  as 
well  as  in  estimating  the  sediment  discharge  at  various  locations  of 
waterway  system.  A  large  number  of  contributions  in  the  area  of  sediment 
transport  in  suspension  has  been  published  in  recent  years  using  the 
numerical  techniques  to  solve  the  sediment  convection-diffusion  equation  in 
a  vertical  plane.  Some  typical  examples  may  be  found  in  Jobson  and  Sayre 
(1970),  Yang  and  Sayre  (1971),  and  Chen  (1971).  Smith  and  O'Connor  (1977) 
and  Kerssens,  et  al.  (1977)  presented  their  findings  at  the  17th  Congress 
of  IAHR.  The  former  paper  described  a  two-dimensional  model  in  a  vertical 
plane  which  gives  the  velocity,  as  well  as  the  concentration  fields,  of  an 
estuarial  type  flow  with  only  good  agreement  with  experimental  data  of  the 
velocity  distribution.  The  latter  paper  succeeded  in  combining  the  quasi¬ 
steady  fluid  flow  equations,  sediment  continuity  equation,  and  the 
convection-diffusion  equation  for  morphological  computations  in  a  vertical 
plane  of  a  very  wide  alluvial  channel.  Some  of  its  basic  assumptions  are 
adopted  from  an  early  work  of  deVries  (1973).  Their  numerical  estimation 
of  sand  bed  deformation  considering  only  sediment  transport  in  suspension 
compares  quite  well  with  their  own  experimental  results.  These  two  models 
are  based  on  finite-difference  schemes.  Leimkuhler  et  al .  (1975),  and 
Ariathurai  et  al.  (1976)  applied  the  FEM  to  obtain  solution  of  the  sediment 
convection-diffusion  equation  in  a  vertical  plane  with  some  success.  Most 
recently,  Alonso  and  Wang  (1980)  presented  the  results  of  a  study  of  local 
scour  and  fill  in  sand  bed  stream.  The  bed  deformation  is  verified  with 
experimental  results.  This  latest  model  is  also  a  two-dimensional  one 
applying  only  to  a  vertical  plane  or  a  case  of  a  very  wide  channel. 
Recently,  a  very  comprehensive  water  and  sediment  routing  model  based  on 
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the  depth  averaging  approach  has  been  developed  by  Simons  et  al.  (1979). 
They  verified  the  applicability  of  an  analytic  model,  supplemented  by 
empirical  relations,  for  simulating  water  and  sediment  movement  in  a  river 
system.  This  model  uses  a  finite-difference  numerical  scheme. 

Both  of  the  previous  approaches  have  been  adopted  in  the  present  study 
to  develop  two  different  finite-element  schemes.  One  uses  a  two- 
dimensional  vertical  domain;  the  other  employs  a  two-dimensional  depth¬ 
averaging  solution.  Detailed  information  on  Mathematical  Formulation, 
Numerical  Modeling  and  Solution,  and  Simulation  Results  are  presented  in 
Parts  2,  3,  and  4,  respectively.  A  complete  computer  program  listing  is 
given  in  the  addendum  2. 

2  FORMULATION  OF  MATHEMATICAL  MODEL 

2.1  EQUATIONS  OF  WATER  MOVEMENT 

In  order  to  make  mathematical  modeling  a  possibility,  many  basic 
assumptions  and  simplifications  are  necessary.  Since  there  is  no  theory 
which  can  include  the  fluid  flow  and  moving  boundary  of  the  sediment 
particles  'simultaneously',  the  flow  in  an  alluvial  channel  is  to  be  broken 
into  two  physical  phenomena  and  studied  in  a  alternating  sequence.  That 
is,  the  hydrodynamic  characteristics  of  a  fluid  flowing  along  a  channel 
with  an  "instantaneously  fixed"  sand  bed  of  a  given  geometry  are  studied 
first,  and  then  the  deformation  of  the  sand  bed  is  calculated  using  local 
sediment  discharge  determined  from  the  hydrodynamic  characteristics  as  well 
as  from  the  sediment  properties.  The  time-dependent  phenomena  of  fluid 
flow  and  bed  deformation  are  simulated  by  carrying  out  these  two  steps  of 
solution  procedure  repeatedly. 

In  developing  the  mathematical  model  to  describe  the  flow  of  water 
along  a  channel  with  "instantaneously  fixed"  sand  bed  boundaries,  the 
conservation  laws  of  mass  and  momentum  for  incompressible,  viscous  fluids 
are  applied.  Written  in  tensor  form,  they  are: 

v.  .  =  0,  i  =  1,  2,  3,  (1) 

X  ,  1 

vi  +  (viv0  i  +  (p  i  "  Tii  -  Fi  3  &  j  =  1,  2,  3,  (2) 

Where  is  the  ith  component  of  the  velocity  vector;  p,  p,  and  v  are 

pressure,  density,  and  kinematic  viscosity  of  the  fluid,  respectively; 
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are  components  of  the  stress  tensor  which  represent  laminar  or 

turbulent  stresses;  F^  is  the  ith  component  of  the  external  (gravitational) 

body  force  per  unit  mass;  v  .  and  v  represent  partial  derivatives  of  the 

» J 

function  v  with  respect  to  the  coordinate  direction  x^  and  Lime,  L, 
respectively.  The  summation  convention  is  adopted  for  repeated  indices 
throughout  this  work.  For  laminar  flow,  the  expression  for  stress  tensor 
is 


T  .  .  =  u(v.  .  +  v.  .) 

ij  i,J  J,i 


(3) 


In  the  case  of  turbulent  flow,  however,  a  different  closure  scheme  is  used 
to  complete  the  mathematical  model.  It  is  beyond  the  scope  of  this  study 
to  review  all  existing  schemes.  Therefore,  only  those  used  in  the  present 
work  will  be  discussed  whenever  they  are  introduced. 

Boundary  conditions  needed  for  the  fluid  flow  simulation  are  that  (1) 
neither  slip  nor  seepage  are  allowed  at  the  channel  boundaries,  where  the 
pressure  is  left  to  be  determined  by  applying  the  governing  equations  at 
the  boundary;  (2)  on  the  free  surface  the  pressure  is  taken  as  constant; 
(3)  the  shear  stresses  acting  on  the  free  surface  are  neglected,  which 
implies  that  the  maximum  velocity  occurs  at  the  free  surface;  (4)  the  flow 
is  considered  to  be  fully  developed  or  uniform  at  the  upstream  end  of  the 
channel,  and  (5)  appropriate  boundary  conditions  are  imposed  at  the 
downstream  end.  Also  assumed  is  that  the  sand  bed  will  be  of  uniform 
roughness  and  fixed  instantaneously.  And,  although  the  sand  bed  has  been 
assumed  flat  to  start  the  simulations,  any  prescribed  bed  form  (not 
necessarily  flat)  can  be  used  as  an  initial  condition  without  difficulty. 
More  about  bed  deformation  will  be  discussed  later.  Now,  the  attention  is 
still  centered  on  the  mathematical  modeling  of  the  open  channel  flow  with 
instantaneously  fixed  boundaries. 

Although  realistic  mathematical  models  of  flows  in  natural  streams 
should  be  both  three-dimensional  and  time-dependent,  the  computer  storage 
capacity  required  to  store  the  information  and  the  computing  time  needed  to 
obtain  converged  solution  is  too  expensive  to  justify  its  use  for  a 
preliminary  analysis  of  basic  characteristics  of  flows  in  an  alluvial 
channel.  Therefore,  a  more  practical  alternative  is  needed.  As  mentioned 
above,  two  viable  alternatives  discussed  in  this  report  are  used  to  treat 
horizontal  shallow  flows,  and  vertical  flows  with  negligible  variations  in 
the  lateraL  direction.  These  two  approaches  are  presented  below. 
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a.  Depth-Averaging  Model:  This  approach  has  been  found  to  he  both 
adequate  and  efficient  for  shallow  water  fiows  in  which  the  variation  of 


the  hydrodynamic  characteristics  is  predominantly  horizontal.  That  is, 
velocity  variations  in  the  horizontal  plane  are  greater  than  variations  in 
the  vertical,  and  the  pressure  distribution  along  the  vertical  axis  may  be 
assumed  to  be  hydrostatic.  In  such  cases,  an  efficient  as  well  as  adequate 
approximation  to  the  three-dimensional  flow  problem  is  to  treat  the 
vertically  averaged  flow  properties  in  a  horizontal  plane.  To  derive  the 
model,  let  us,  first,  represent  the1  bee!  and  free  surface  geometries  by 
f(x,y,t)  and  f)(x,y,t)  respectively  as  shown  in  Tig.  1. 

Although  both  the  channel  bed  and  tree  surface  of  the  flow  are,  in 
general,  deformable  in  time,  for  the  convenience  of  numerical  solution, 
they  are  assumed  frozen  instantaneously  during  the  solution  of  flow 
properties.  But,  when  the  flowfield  is  solved,  the  elevations  of  both 
channel  bed  t(x,y,t)  and  free  surface  of  the  flow  r)(x,y,tj  are  replaced  by 
their  new  values  before  the  flowfield  is  solved  again  for  the  next  time 
step.  The  method  to  obtain  new  values  of  t.  and  rj  will  be  given  later  when 
the  technique  for  estimating  sediment  discharges  is  described. 

The  boundary  conditions  described  in  the  previous  section  can  then  be 
translated  into  the  mathematical  equations.  On  the  instantaneously  frozen 
channel  bed,  z  =  f(x,y),  they  are: 

v.  =  0,  i  =  1  ,  2  and  '1  (A) 
and  at  the  free  surface,  z  =  !"|(x,y,t),  they  are: 


p  =  p 

o  *  v°  r)  j  =  Vv 
t  . ,,  -  t  n  .  -  i's , 

i  1  i.j  ,J  ' 
Where  p  is  the  atmospheric 


j  =  1  and  2 
i  and  j  =  1  and  2 


(5) 

(6) 

(7) 


pressure  and  t‘.  is  the  i  th  component  of  the 
A  similat  expression  may  he  written  for  the 
bed  shear  stress,  t”.  Integrating  Kqs .  1  and  2  and  using  U  and  V  to 

if  Vj  and  v^  respectively,  one 


surface  (wind)  shear  stress. 
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where  g  is  the  acceleration  of  gravity, 
acting  on  vertical  planes.  During  the 


and  T . .  are 
ij 

integration , 


effective  stresses 
the  Leibnitz  rule, 


e-g- 


3  n  n  cif 

f  f(x,  y,  x,  t)  dz  =  /  gi  dz 
+  f(x,  y,  n,  t)  -  f (x,  y,  £,  t)  |j> 


(12) 


and  the  boundary  conditions  (6)  and  (7)  have  been  applied. 

The  bed  shear  stresses  are  assumed  to  have  the  same  magnitude  as 
those  in  steady  uniform  flow,  and  their  directions  to  be  the  same  as 
those  of  the  depth-averaged  velocity  components.  Their  mathematical 
representations  are 

(i^  ;  tb  =  pf(U2+V2)*  (U;V)  (13) 

A  y 

where  f  is  a  dimensionless  friction  factor  defined  in  terms  of  either 
the  Chezy  coefficient,  C^.,  or  the  Manning’s  roughness  coefficient,  n, 
as  given  below: 


f  =  g  C{*  or  f  =  g  n  (1.49)  h  '  (14) 

To  enable  the  model  to  simulate  eddies  and  circulations  in  the  depth- 
averaged  two-dimensional  flow,  the  effective-stress  terms  have  been 
taken  into  account.  By  analogy  with  the  eddy-viscosity  approach  used 
in  some  turbulence  closure  schemes,  an  eddy-viscosity,  £,  is 
introduced,  so  that  the  effective-stress  terms  in  Eqs .  8  and  9  are 

replaced  by 


92U . 
9y^ 


(15) 
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and , 


respectively.  As  a  first  approximation,  the  coefficient  of  eddy  viscosity 
can  be  related  to  the  mean  flow  properties  by 

e  =  Eh  (U2  +  V2)*5  (17) 

where  E  is  an  adjustable  eddy- viscos i ty  parameter.  Thus,  the  derivation  of 
the  depth-averaging  model  has  been  completed. 

b.  Two-Dimensional  Vertical  Flow  Model:  In  this  case  lateral  changes  in 
flow  and  bed  material  properties  are  neglected.  In  tracking  the  bed- 
profile  evolution,  advantage  is  taken  of  the  fact  that  the  bed  deforms  at  a 
rate  much  slower  than  the  rate  of  change  of  free-surface  transients.  The 
water  discharge  hydrograph  is  replaced  by  a  piecewise  continuous 
distribution  with  time  increments  smaller  than  the  time  scale  of  the  bed 
transients.  During  these  time  increments  the  water  surface  profile  is 
computed  by  assuming  a  one-dimensional  spatially  varied  steady  flow.  This 
permits  the  uncoupling  of  the  bed  profile  calculations  from  the  water 
routing  scheme.  Thus,  the  flow  governing  equations  are  the  well  known 
spatially  varied  flow  equations 

U  U,x  +  8  h  x  +  8  Zb  x  +  8  11  fu!/Cf  h  "  0  O8) 

and 

u  h  =  q  (19) 

where  u  is  the  depth-averaged  velocity,  h  is  the  flow  depth,  z^  is  the  bed 
elevation,  x  is  the  streamwise  distance,  g  is  the  acceleration  of  gravity, 
Cj.  is  the  Chezy  friction  coefficient,  and  q  is  the  water  discharge  per  unit 
width.  Within  the  context  of  the  one-dimensional  flow  assumption,  the  flow 
geometry  is  restricted  to  situations  where  no  separation  occurs  and  where 
nearly  parallel  flow  conditions  exist.  Under  these  conditions  the  vertical 
velocity  distribution  at  any  section  is  assumed  to  follow  the  logarithmic 
1  aw 

ii.,.  10  (z-z,  ) 

u  =  ~  In  -  k  (20) 

s 


-eh  ( 


a2v 

3x2 


a2v, 


(16) 


Kit 


in  which  u.,.  is  the  local  bed  shear  velocity;  k  is  the  von  Karman  constant, 
k  is  the  equivalent  grain  roughness,  and  z  denotes  the  vertical  position. 
Integrating  Eq.  20  over  the  flow  depth  gives 

uv  =  K  il  (  In (10  h/ks )  -ll'1  (21) 
which  permits  the  estimatiot.  of  u.,_  from  known  local  flow  parameters. 

2.2  EQUATIONS  OF  SEDIMENT  MOVEMENT 

As  discussed  in  the  previous  chapter,  the  sediment  transport 
phenomenon  associated  with  the  flow  in  an  alluvial  channel  is, 
unfortunately,  too  complex  to  be  described  completely  by  analytical 
techniques,  because  the  theories  lor  describing  the  bed  load  of  sediment  as 
well  as  the  boundary  condition  of  sediment  concentration  on  the  bed  surface 
have  not  yet  been  developed.  The  several  practical  sediment  transport 
models  that  have  been  adopted  most  often  can  be  classified  into  three 
categories : 

1.  By  assuming  that  the  sediment  is  being  transported  predominantly  in 
suspension,  the  bed  load  is  included  in  the  suspended  load  and  the  sediment 
diffusion-convection  equation  can  be  solved  for  the  sediment  concentration 
distribution  by  analytical  or  numerical  schemes.  Then,  the  sediment 
discharge  is  computed  by  integrating  the  mass  flux  over  a  cross-section 
area  normal  to  the  flow  direction.  The  drawbacks  of  this  approach 
are:  first,  the  boundary  condition  for  sediment  concentration  on  the 
surface  of  the  channel  bed  is  still  not  understood  completely.  Secondly, 
the  estimation  of  sediment  discharge  may  not  have  sufficient  accuracy  for 
the  cases  in  which  the  bed  load  can  no  longer  be  neglected. 

2.  The  sediment  discharge  is  estimated  by  combining  the  suspended 
discharge  and  the  bed  discharge.  The  suspended  discharge  is  solved 
analytically  using  the  approach  similar  to  that  described  above  in  the 
category  i.  The  bed  discharge  is,  however,  estimated  using  one  of  the 
empirical  functions  which  have  been  developed  in  recent  years.  A  good 
review  concerning  the  capabilities  and  limitations  on  those  bed  load 
form  ilas  may  be  found  in  Graf  (1971).  The  major  drawbacks  of  these  kind  of 
approaches  are  that,  first,  they  are  more  involved  than  the  first  category 
approach,  and  second,  it  is  difficult  to  find  an  accurate  formula  for 
determining  the  concentration  boundary  condition  near  the  surface  of  the 
cl i a r. ne 1  bed.  This  latter  dr'whack  has  to  be  resolved,  before  this  approach 

to  he  widely  adopted. 


3.  Without  distinguishing  the  suspended  and  bed  load,  one  may  select  one 
appropriate  empirical  formula  to  compute  the  total  sediment  transport  load. 
The  obvious  advantage  is  the  simplicity  of  this  approach  when  it  is  used  to 
simulate  a  particular  section  of  a  particular  alluvial  channel.  However,  a 
different  formula  may  have  to  be  adopted  for  a  different  channel  or  even  a 
different  section  of  the  same  channel.  This  drawback  can  sometimes  be 
alleviated  by  calibrating  an  adjustable  parameter  existing  in  the  empirical 
function  chosen  for  a  particular  channel,  or  even  calibrating  it  from 
section  to  section  of  a  channel,  so  that  an  accurate  sediment  transport 
model  can  be  established. 

Although  all  three  different  approaches  have  been  used  in  the  course 
of  this  research,  the  mathematical  formulations  of  the  first  and  the  third 
approaches  are  described  below,  because  they  have  produced  some  good 
results . 

a.  Depth-Averaging  Model:  In  simulating  depth-averaged  flows 
predominantly  in  a  horizontal  plane  the  total  sediment  transport  load  is 
represented  by  an  empirical  formula.  For  example,  a  model  based  on  a  power 
function  of  the  mean  velocity  may  be  used  conveniently  to  estimate  the 
total  load  with  its  coefficient  and  exponent  being  calibrated  for 
individual  cases.  This  model  has  been  successfully  adopted  by  other 
researchers  (Simons  et  al.,  1979,  DeVries,  1973).  Therefore,  it  is  applied 
to  perform  a  preliminary  study  on  the  deformation  of  sand  bed  due  to 
sediment  transport  by  the  fluid. 

The  sediment  transport  function  selected  for  the  present  study  is  of 
the  form 


8i 


c  p  V1" 
g  i 


(22) 


where  g^  is  the  bed  material  transport  function  in  the  direction  of  depth- 
averaged  velocity  component  ,  c^  and  m  are  the  coefficient  and  exponent 


respectively.  For 


particular  channel  geometry  and  bed  material 


characteristics,  c  and  m  are  estimated  or  calibrated  using  experimental 
8 

information.  After  the  depth-averaged  flow  field  is  solved,  the  sediment 
discharge  can  be  easily  computed. 

Then,  the  bed  material  continuity  equation 


q  +  q  +  (1  -  \)  z,  =  0 
^s,x  s  ,y  b,t 


(23) 


is  used  to  calculate  a  new  bed  elevation  z^.  And,  the  complete  physical 
system  has,  thus,  been  represented  by  a  mathematical  model.  In  this 
equation  \  denotes  the  porosity  of  the  bed  material. 
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b.  Two-Dimensional  Vertical  Flow  Model:  The  preceding  profile 
computations  are  based  on  the  assumption  that  the  local  flow  transport 
capacity  does  not  deviate  significantly  from  equilibrium  conditions,  and  no 
distinction  is  made  between  bed  and  suspended  transport  mode.  Such  an 
approach  is  justified  if  either  bed-load  transport  is  predominant,  or  the 
sediment  is  being  routed  over  long  reaches.  Otherwise,  the  assumption  of 
near  equilibrium  capacity  is  not  valid  because  the  bed  load  reacts 
immediately  to  changes  in  local  flow  conditions  while  the  suspended  load 
tends  to  react  more  slowly.  This  slow  adaptation  of  the  suspended  load  can 
significantly  influence  the  bed  profile  evolution  over  relatively  short 
reaches  (Kerssens  et  al.,  1977);  Freds*$e,  1978).  For  this  reason,  the 
present  analysis  includes  the  equation  governing  the  dispersion  of 
suspended  sediment,  in  addition  to  the  bed-material  continuity  equation. 
These  equations  are  simplified  by  assuming  that: 

i.  Longitudinal  dispersion  can  be  neglected  in  relation  to  the  vertical 
dispersion; 

ii.  The  sediment  settling  velocity,  v  ,  is  invariant; 

iii.  Vertical  convection  is  negligible  in  a  nearly  horizontal  fiow; 

iv.  The  time  rate  of  change  of  the  sediment  concentration  is  not 
significant  (Mahmood,  1975); 

v.  The  bed  material  is  fairly  uniform  in  size,  and  can  be  represented  by 
an  effective  particle  diameter. 

The  sediment  dispersion  equation  thus  yields 


uc  -vc  —  ( £  c  ) 

,  x  s  ,  z  z  ,  z  ,  z 


(24) 


where  c  is  the  point  volumetric  concentration,  and  e  is  the  vertical 
sediment  transfer  coefficient.  By  averaging  Eq .  24  over  the  flow  depth  the 
following  bed-material  continuity  equation  is  obtained 


qs,x  +  (1"X)  Zb,t  =  ° 


(25) 


in  which  qg  is  the  volumetric  discharge  of  bed  material  per  unit  width,  K 
's  the  bed  porosity,  and  t  is  time. 

Since  Eqs .  18,  19,  20,  24,  and  25  describe  an  evolutionary  process, 
appropriate  initial  and  boundary  conditions  need  to  be  specified.  They  are 
the  initial  bed  profile,  the  upstream  bed  elevation  and  inflow  of  water  and 
sediment,  the  downstream  flow  stage  and  streamwise  concentration  gradient, 
and  vanishing  sediment  flux  across  the  free  surface.  These  conditions  are 
represented  by 
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zb(x;0)  =  Zq ( x )  ,  0  $  x  <1M  t  =  0  (26) 
q(0;t)  =  g(t)  ,  x  =  0,  t  ^  0  (27) 
qs(0;t)  =  r(t)  ,  x  =  0,  0,  (28) 
c(0,z;t)  =  s ( z ; t ) ,  x  =  0,  0§z-zb$h,  t^O  (29) 


zfa  +  h  =  H  (t)  ,  x  =  L,  t  a  0, 


(30) 


C,x  ~  KV.z3*2  +  Vs  r,z1/u>  x  "  L>  0  -  z  '  zb  --  h,  t  *  0,  (31) 
vgc  +  tz  c  z  =  0  ,  0  i  x  $  L,  z  =  zb  +  h,  t  £  0.  (32) 


In  these  equations,  zQ,  g,  r,  s,  and  H  are  given  continuous  functions,  H 

represents  the  downstream  stage,  and  L  is  the  length  of  the  reach.  The 

sediment  concentration  over  the  upstream  boundary,  condition  29,  is 

determined  by  solving  Eq.  24  for  the  equilibrium  case  (c  =  0)  and  by 

> x 

specifying  an  appropriate  distribution  for  t  .  From  regression  analysis  of 
point  measurements  of  the  sediment  transfer  coefficient  conducted  by 
Coleman  (1970),  Kerssens  et  al.  (1977)  obtained  the  following  expression 
for  the  transfer  coefficient. 

e  =  e  [4  £  (1-t)]6  ( 33) 

where 

£  =  (z-zb)/h, 

£  =  u,.h  (0.13  +  0.2  (v  /  u... ) 2  •  1  2  ] 

max  -  s  •' 

6  =  1,  0  <  C  ^  0-3, 

6=0,0. 5  <£Si. 

Using  this  equation  the  following  equilibrium  concentration  distribution 
results 

c  =  hu*c  [t  /(1-t  )]P  x  l(l-t,)/t.]p6exp[2p(6-l)(2t-l)],  (34) 

in  which 


P  =  h  v  /4e 

s  max 

and  c  is  the  concentration  at  a  point  C,  >0  near  the  bed.  An 
a  a 

the  reference  value  cg  is  obtained  by  assuming  that  most  of  the 
carried  in  suspension,  and  then  matching  the  sediment  flux 
upstream  section  to  a  suitable  total  load  formula.  That  is, 


estimate  of 
sediment  is 
through  the 
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(35) 


1 

q  =  /  u  c’df 
S 

Kerssens  et  al.  (1977)  found  the  following  modified  version  of  the  formula 
developed  by  Engelund  and  Hansen  (1967)  to  agree  well  with  their  measured 
loads 

qs  =  0.035  uS  /  [d5Q  (S- 1 ) 2  g  ^  C|)  ,  (36) 

where  q  is  the  sediment  volumetric  discharge  in  m3/sec.m,  d^„  is  the  bed 
s  50 

material  size  in  meters,  and  S  is  the  sediment  specific  gravity. 
Substituting  Eqs .  20,  34  and  36  into  Eq.  35  yields 

— ^ —  1(1  -  1  )/£  JP  = 

hu,.d50(S-l)2g'5C3  3  -a 

c  J£n(30£/k  )[(l-£)/UP6  exp  (2p(6-l )  (2£-l)  Ut  (37) 

d  ^  o 

''a 

from  which  c  can  be  determined.  The  integral  on  the  right  side  is 

a 

evaluated  using  Gaussian  quadrature. 
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NUMERICAL  SOLUTIONS 


3 

3.1  DEPTH-AVERAGING  MODEL 

A  typical  finite  element  system  used  to  discretize  a  continuous 
horizontal  domain  is  shown  in  Fig.  2. 

Isoparametric  interpolation  functions  are  used  to  relate  the 
properties  at  any  point  within  an  element  to  those  at  nodal  points  of  the 
same  element.  For  linear  cases,  they  are  of  the  form: 


=  ®n  (U„;  Vn;  h,) 


where: 


<*n  =  k  (i+£nO  (i+nnn),  n 


=  1 ,  2 ,  3 ,  and  4 


The  transformation  between  the  global  (x,y)  and  the  local  (£,i))  coordinate 
systems  is  specified  by 

(X;  Y)  =  «n  (Xn,  Yn).  (40) 

To  reduce  the  degree  of  nonlinearity,  the  hydraulic  depth,  h,  in  the 
momentum  equations  at  all  nodes  are  estimated  initially  and  corrected 
during  the  iterative  solution  procedure.  Due  to  the  fact  that  the  bed 
deforms  at  a  rate  much  slower  than  the  rate  of  change  in  hydrodynamic 
properties,  the  hydrodynamic  equations  and  the  bed  form  model  are 
decoupled.  And,  since  the  bed  form  can  be  considered  being  instantaneously 
fixed  during  the  solution  of  hydrodynamic  equations,  the  assumption  of 
quasi-steady  flow  should  give  satisfactory  solution.  Therefore,  the  fluid 
flow  equations  are  simplified  into  the  following  form: 


U  +  y  9U 

ax  dY 


,a2u  ^  a2u.  .  gw*  . .  ,  an 

^  gy £  ^  f>  av 


!>  ♦  cS  u  *  * 


..  av  .  ..  av  c  ,a2v  a2 
u  ax  v  by  ‘  (axTf  3Y 


,3 2v  A  a2v,  .  gw*  „  .  an  _ 
(aX*  3Y1-1  c|h*  V  8  3Y  ° 


.9u  ,av  .ah  ah_ 
h  3X  h  3Y  u  ax  v  3Y  0 


(41) 

(42) 

(43) 


where  W*  =  (U*2  +  V*2)*6 

Applying  the  Method  of  Weighted  Residual  (MWR)  after  Galerkin,  the 
finite  element  equations  for  each  and  every  element  are  formulated. 
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0 


where 


A  0  U.U  +  B  V.U  +  C  U  +  D  U  +  E  h  -  H  = 
n£m  £  m  n£m  £  m  rim  m  11m  m  nm  m  n 


A  .  U  .V  +  B.V.V  +  C  V  +  D  V  +  F  h  -  G  =  0 
n£m  £  m  n£m  £m  nmm  nmm  nm  m  n 


A  .  h  U  +  B  .  h.V  +  A  .  U.h  +  B  .  V.h  =  0 
n£m  £  m  n£m  £  m  n£m  £  m  n£m  £  m 


A  .  =  fa  ana  dA 

n£m  .  n  £  m,x 
A  ’ 


B  .  =  fa  ana  dA 

n£m  J  n  £  m,y 


C  =  fe(oi  a  +  a  a  )dA 
nm  J  n,x  m,x  n,y  m,y 


:W*  r 

■x  J  a  a  dA 


nm  C‘h*  n  m 
f  A 


E  =  gja  a  dA 

nm  ^  r>  m  n 

A 


F  =  g/a  a  dA 
nm  ^  n  m,y 


at 


H  =  Jea  (U  n  +  U  n  )d£  -  e  ^  Ja  dA 
n  J£  nv  ,x  x  ,y  y'  s  dx  ^  n 

— 

G  =  Jea  (V  n  +  V  n  )d£  -  g  r"  J  a  dA 
n  £  n  ,x  x  ,y  y  6  3y  ^  n 


(44) 

(45) 

(46) 


All  sets  of  finite  element  equations,  one  set  for  each  element,  are 
assembled  into  the  following  global  set.  The  actual  assembling  process  is 
carried  out  by  a  computer  program  simultaneously  with  the  evaluation  of  the 
coefficients  of  local  element  equations. 


. u.u, 

+ 

B.  V.U, 

+ 

C  . ,  U,  +  E  .  . 

h  -  H . 

(47) 

ijk  ,1  k 

ijk  j  k 

lk  k  i j 

J  i 

•  u.v. 

+ 

B.  V.V, 

+ 

C..  V.  +  F.  . 

h.  -  G. 

(48) 

ijk  J  k 

ijk  J  k 

lk  k  ij 

J  i 

.  ..  h.U, 

+ 

B.  h.V, 

+ 

A.  ..U.h,  + 

B  V.h. 

(49) 

ijk  J  k 

ijk  .]  k 

1  jk  j  k 

ijk  j  k 

This  global  set  is  solved  iteratively  by  the  Newt on-Raphson  technique 
to  yieLd  the  flow  field  properties. 

Then  the  sediment  transport  function,  F.q .  22,  is  used  to  estimate  the 
bed  material  discharge  at  eacn  node.  The  sediment  continuity  equation  is 
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discretized  by  the  FEM,  using  the  same  element  system,  and  interpolation 
functions.  The  resulting  global  equations  are: 


L .  .  AZ  .  =  K .  .  (50) 
ij  j  1 

Substituting  the  sediment  discharge  results  into  this  equation  the  new  bed 
elevations  at  every  node  are  obtained.  Thus,  the  discretization  of  the 
mathematical  model  of  the  sediment  transport  phenomenon  has  been  completed. 
3.2  TWO-DIMENSIONAL  VERTICAL  FLOW  MODEL 


The  Galerkin  integral  approximation  within  the  Method  of  Weighted 
Residuals  is  used  to  derive  the  finite  element  equations.  The  first  step 
of  the  solution  involves  dividng  the  continuous  solution  domain,  0,  into  a 
set  of  nonoverlapping  subdomains,  or  finite  elements,  such  that 

E 

fi  U  T  =  U  Q 
,  e 


where  f  is  the  boundary  of  the  solution  domain,  and  E  denotes  the  total 
number  of  elements.  The  one-dimensional  domain  of  Eq .  18  is  discretized  as 
shown  in  Fig.  3(a).  In  Fig.  3(b)  the  two-dimensional  domain  of  Eq.  24  is 
divided  into  a  series  of  isoparametric  elements.  The  finite  element  grid 
is  designed  so  as  to  provide  greater  resolution  near  the  stream  bed  where 
high  velocity  and  concentration  gradients  exist.  Within  each  subdomain  the 
unknown  variables  are  interpolated  in  terms  of  their  nodal  values  as 

fg(x;t)  =  4>n(x)  Fn(t)  ,  x  e  0p>  (51) 

where  f  stands  for  the  interpolated  values  of  all  variables,  F  denotes 
e  r  n 

values  at  a  node  n,  and  0^  represents  the  interpolating  functions.  This 
study  uses  linear  isoparametric  functions,  which  have  the  advantage  of 
satisfying  the  basic  convergence  criteria  of  completeness  and  continuity, 
and  serve  at  the  same  time  as  coordinate  transformation  equations.  For 
instance,  for  the  two-dimensional  elements 

<t>n  =  k  (i+tnt)0+nn  n)  (52) 

where  the  isoparametric  coordinates  are  related  to  the  global  coordinates 
such  that  the  corner  point  (x.  ,  z  ,  i  =  1,  2,  3,  4)  of  a  quadrilateral 

element  are  transformed  into  the  four  points  (1,1),  (1,-1),  (-1,-1),  (-1,1) 
in  the  (t,q)  space.  The  approximate  solutions  are  of  the  form 

E 

f(x;t)  =  U  f  =0  1-  ,  x'r.Q,  (5.3) 

,  e  n  it 
e=l 
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tch  of  finite  element  grid  at  t  =  0 


in  which  4>^  are  Linear  functions  of  the  spatial  domain.  When  the  solutions 
?  are  substituted  into  the  Eqs .  18,  19,  and  24  approximation  errors,  or 
residuals,  are  introduced.  The  expansion  coefficients  F  are  determined  by 
rendering  these  residuals  orthogonal  to  the  respective  interpolation 
functions.  This  yields 


J  [g-q2/h3)h  +  g(z  +  q2/C2h  !)  ]<{>  dQ  =  0 

n  >  *  u  > x  1  n  e 


.17 


Q 

e 


(v 


£  ) 
Z  ,Z 


C  -£  C 
,z  z  ,zz 


4>  dQ  =  0 
n  e 


(54) 

(55) 


The  integrand  in  Eq.  54  results  from  combining  equations  18  and  19.  After 
substituting  the  interpolation  functions,  applying  the  Green-Gauss  theorem 
to  the  last  term  in  Eq.  55,  or  integrating  by  parts,  the  preceding 
equations  yield 


L  H 
nm  m 


Mn> 


P  0  U0C  +  Q  C  =  R  , 
n£m  £  m  nm  m  n 


(56) 

(57) 


where  £,  m,  n  =  1,  2,...,  n^  (n^  =  total  number  of  element  nodes),  and 


L  =  /  (g-q2/h3)  <J>  <p  dfi  , 

nm  q  ^  n  m,x  e’ 


Mn  =  g(zb,x  +  q"/Cf  f‘3)  W 
e 

p  „  =  JJ  0  <t>  „ct  da  , 

n£m  ^  n  £  m,x  e 

e 

0  =  JJ  e  (<{>  4>  -  £  (J)  <!)  i  da 

‘nm  q  z  n,z  m,z  z,z  n  mmz  e 

e 


(58a) 

(58b) 

(58c) 

(58d) 


R  = 
ri 

ft  <  c  <t>  ) 

jl  z  ,  z  ,  z  n ' 

e 

n  dl  . 
z  e 

(58e) 

Here  n  is 

z 

the 

z-component  of 

the  unit 

vector  normal  to 

the 

element 

boundary  T  , 

and 

h,  c  represent 

either  ini 

tial  or  iterated 

values.  The 

coe f f i c i en ts 

(58) 

are  evaluated 

in  the  (f, , 

r| )  space.  Since 

all 

elements 

are  identical  in  this  coordinate  system,  once  the  formulas  for  those 
coefficients  are  developed  for  one  element  the  same  formula  is  used  for  all 
the  other  elements.  Eq .  (58c)  yields,  for  instance, 


> 
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(59) 


P  „  =  IS  4>  0 df) 

nJ£m  n  £  m,x  e 

e 
+  1 


=  {{  Vi  •♦.,cJn*  VnJi2)  |J|  dt'dn 


in  which  [Jl  is  the  determinant  of  the  transformation  Jacobian,  and  J..  are 

ij 

elements  of  the  inverse  Jacobian.  In  all  cases  the  integrals  are  evaluated 
using  a  Gaussian  quadrature  subroutine. 

Once  evaluated,  the  coefficients  of  equations  56  and  57  are  assembled 
into  the  corresponding  global  matrices.  After  introducing  the 
Dirichlet-type  boundary  conditions,  the  global  equations  reduce  to 


iu  h  ‘  \  •  (60> 

Tki  C£  =  Rk  •  (61) 

where  T^£  =  P^£  +  Qj^,  an<^  ^ >  £  =  1>  2,...,  N.  Here  N  is  the  total 

number  of  node  points  containing  unknown  values.  The  nonlinear  equations 

60  and  61  are  solved  iteratively  using  a  subroutine  developed  for  banded 

matrices.  After  the  global  nodal  values  are  determined,  Eq.  19  is  used 

to  obtain  the  corresponding  velocities  U^. 

For  convenience,  the  sediment  continuity  Equation  25  is  discretized 

using  the  following  explicit  finite-difference  scheme: 


i»j  +  l 


i  .J 


q  i  +  1’j  -  q  iM*j 
At  4s  Hs 


O-A)  x.+1  -  x._1 


(62) 


where  x..,  and  x.  ,  are  defined  in  Figure  1(a),  and  At  =  t.  ,  -  t.  is  the 

l+l  i-l  6  j+1  J 

time  increment.  Equations  60  through  62  constitute  the  numerical  algorithm 
used  in  the  present  case. 
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4  PRESENTATION  AND  EXAMINATION  OF  RESULTS 

4.1  FLOW  AROUND  A  SPUR  DIKE  IN  A  STRAIGHT  CHANNEL  WITH  A  SAND  BED 

A  thin  spur  dike  is  placed  in  a  straight  alluvial  channel  with 
dimensions  shown  in  Fig.  4.  The  variation  of  the  flow  field  hydrodynamic 
behavior  due  to  bed  form  changes  in  time  is  simulated  using  the  depth¬ 
averaging  model  developed  in  the  present  study.  The  boundary  as  well  as 
initial  conditions  used  in  this  case  are  no  slip  on  the  solid  walls,  almost 
uniform  flow  at  the  entrance,  parallel  flow  are  left  to  be  determined  by 
the  governing  differential  equations  at  the  exit,  no  shear  on  the  free 
surface  of  the  water,  bed  material  particles  having  uniform  properties,  and 
the  bed  being  flat  when  the  flow  simulation  is  started.  The  computer 
simulated  bed  form  and  velocity  distribution  in  the  flow  field  after  5,  10, 
and  15  hours  are  shown  in  Figs.  5  through  10  for  an  average  velocity  = 
0.110  ft/sec;  water  depth  =  0.38  ft;  channel  width  =  6  ft;  dike  opening  =  4 
ft;  and  the  time  step  for  sediment  routing  was  t  =  3000  sec.  The  velocity 
of  flow  accelerates  locally  increasing  the  carrying  capacity,  which  can  be 
seen  in  the  area  around  the  nose  of  dike  in  Figs.  7  and  10.  When  the  flow 
decelerates  downstream  from  the  dike  and  loses  carrying  capacity,  the 
deposition  of  sediment  occurred.  At  the  downstream  end,  the  flow  recovers 
equilibrium,  and  the  carrying  capacity  is  equivalent  to  the  amount  of 
sediment  deposition.  Therefore,  the  bed  form  essentially  does  not  change. 
This  test  is  a  simulation  of  the  bed  scour  around  a  spur  dike  measured  by 
Zaghloul  and  McCorquodale  (1975).  They  reported  the  geometry  of  the  scour- 
hole  developed  around  a  dike  mounted  perpendicular  to  the  flow  direction  in 
a  laboratory  flume  with  a  movable  bed.  The  maximum  scour  depth  was 
observed  right  at  the  nose  of  the  spur  dike.  The  scour  hole  upstream  of 
the  dike  was  conical  in  shape,  whereas  downstream  was  elongated  and  had  a 
shai lower  slope.  No  sediment  deposition  was  reported.  However,  similar 
measurements  conducted  by  Ahmad  (1953)  and  by  Garde  et  al.  (1961)  have 
shown  that  a  deposition  bar  is  formed  adjacent  to  the  spur  dike  on  the 
downstream  .-id*"1.  The  simulated  depth-average  flow  pattern  and  bed  geometry 
obtained  after  15  (real  time)  hours  exhibit  qualitative  agreement  with  the 
shapes  of  the  measured  hole.  However,  the  locations  of  the  deepest  point 
an  I  deposition  bar  are  incorrect.  Che  reason  is  twofold:  (i)  the  model 
d  .i  not  reproduce  the  flow  separation  created  by  the  stagnation  region 
m  'tr-am  of  the  dike.  The  result  is  an  excessive  flow  component  parallel 
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Fig.  4  Dimensions  of  Straight  Channel  with  Spur 
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to  the  dike  that  forces  the  region  of  maximum  velocity  away  from  the  dike- 
nose;  (ii)  the  low  velocities  predicted  along  the  shear  layer  emanating' 
from  the  dike  induce  a  backflow  circulation  that  is  too  weak  to  transport 
bed-load  material  from  the  primary  stream  to  the  downstream  side  of  the 
dike.  The  results  point  out  a  deficiency  in  the  model  when  simulating 
situations  dominated  by  regions  of  flow  separation.  Further  research  is 
needed  to  correct  this  deficiency.  The  CPU  time  required  for  this 
simulation  was  ten  minutes  per  each  hour  of  real  time. 

4.2  LOCAL  SCOUR  AND  FILL  OF  A  TRENCH  IN  A  SAND  BED  STREAM 

The  natural  backfilling  of  a  trench  dredged  across  an  alluvial 
streambed  is  also  studied  using  the  two-dimensional  model  in  a  vertical 
plane.  This  simplification  is  valid  provided  that  both  the  stream  and  the 
trench  are  wide  enough  to  neglect  the  variation  of  properties  in  the 
lateral  direction  of  the  flow.  The  computer  program  developed  for  this 
simulation  is  basically  similar  to  the  depth-averaged  model  in  a  horizontal 
plane,  so  that  it  is  not  necessary  to  present  it  here. 

The  schematics  of  the  finite  element  system  and  initial  bed  form  and 
flow  conditions  are  shown  in  Fig.  3.  The  dynamic  behavior  of  the  finite 
element  system  is  clearly  seen  from  Fig.  11,  which  represents  the  condition 
of  the  flow  and  deformed  sand  bed  as  well  as  free  surface  at  a  later  time 
(t  >  0).  The  coordinate,  or  position,  of  each  node  is  adjusted  at  each 
time  step  whenever  the  bed  and/or  surface  elevation  change  during  the  time 
interval.  By  doing  this  the  computing  scheme  is  more  stable  and  efficient 
than  the  one  with  fixed  interior  nodes  and  moving  boundary  nodes. 

The  computer  simulated  bed  and  free  surface  profiles  at  7  and  14  hours 
after  the  numerical  experiment  is  turned  on  are  shown  in  Fig.  12.  It  is 
very  obvious  that  the  computer  simulation  can  indeed  produce  reasonable 
results,  which  are  verified  by  data  obtained  from  physical  experiment. 
This  case  was  studied  for  the  purpose  of  verification  with  a  set  of  data 
collected  at  the  de  Voorst  Hydraulic  Laboratory,  The  Netherlands,  by 
Kerssens,  et  al.  (1977).  In  that  experiment,  a  trench  having  the  geometry 
shown  in  Fig.  3  was  formed  in  a  0.127  mm  sand  bed  in  a  laboratory  flume. 
Water  was  circulated  over  the  bed  at  a  steady  rate  of  0.193  m’/sec-m,  while 
the  downstream  stage  was  maintained  at  a  constant  level.  Sediment  was 
supplied  at  a  constant  rate  of  0.041  kg/sec-m,  and  most  of  the  load  was 
carried  in  suspension.  Due  to  hydraulic  sorting  the  d^  of  the  suspended 
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bed  material  approached  a  value  of  0.110  mm  with  a  settling  velocity  of 
O.Olm/sec.  During  the  run  ripples  formed  on  the  bed  and  the  Nikuradse  sand 
roughness  was  estimated  at  0.025  m.  The  Chezy  friction  coefficient, 
corrected  for  side-wall  effects,  was  about  36  m2/sec  on  the  reach  upstream 
of  the  trench. 

The  bed  profile  was  measured  after  7  and  ’4  hours  of  continuous  flow. 
These  profiles  are  shown  in  Figure  15  along  with  the  simulated  free-surface 
and  bed  elevations.  A  reference  concentration  level  equal  to  0.015  m,  or 
slightly  over  half  the  ripple  height,  was  used  in  the  simulation.  The 
comparison  between  the  observed  and  simulated  bed  profiles  is  satisfactory. 
The  largest  shape  discrepancy  occurs  on  the  converging  sideslope  of  the 
trench.  This  may  have  been  caused  by  the  failure  of  the  assumed  velocity 
profiles  to  account  for  the  increase  in  velocity  gradients  near  the  bed 
induced  by  the  flow  acceleration.  The  water-surface  profiles  correctly 
reproduce  the  surface  rise  over  the  trench.  As  the  trench  is  filled  up  the 
water  crest  diminishes  and  moves  downstream  with  the  trench,  as  expected. 

Fig.  13  depicts  the  distribution  of  suspended  sediment  throughout  the 
flow  domain  after  7  hours.  The  contours  of  constant  concentration  were 
drawn  from  computed  results.  This  plot  illustrates  the  deposition  of 
sediment  occurring  on  the  downstream  side  of  the  trench  as  the  stream 
decelerates  and  loses  carrying  capacity.  As  the  flow  regains  velocity, 
sediment  is  entrained  along  the  upstream  side  and  diffused  upwards.  The 
suspended  material  leaving  the  trench  area  gradually  approaches  a  new 
equilibrium  distribution.  The  CPU  time  was  approximately  four  minutes  per 
hour  of  real  time.  The  computer  memory  required  by  the  codes  used  in  the 
tests  presented  above  was  100,000  words. 
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g.  12.  Evolution  of  free-surface  and  bed  profiles. 


5  CONCLUSIONS  AND  RECOMMENDATIONS 

5.1  CONCLUSIONS 

1.  A  two-dimensional  numerical  model  has  been  developed  to  predict  water 
and  sediment  movement  and  water  surface  and  bed-elevation  changes  in 
channel  reaches  with  complex  boundary  geometry. 

2.  The  model  is  based  on  the  conservation  laws  of  water  and  sediment 

continuity,  and  momentum  equations.  The  water  continuity  and  momentum 
equations  are  solved  first.  The  predicted  flow  variables  are  then 

introduced  into  transport  formulas  and  the  sediment  continuity 
equation  to  compute  sediment  load  rates  and  changes  in  bed  elevation, 
respectively.  The  equations  of  motion  are  solved  using  a  finite- 

e 1 emen  t  s  c heme . 

3.  The  model  has  been  validated  by  simulating  laboratory  data  obtained 

from  a  trench  scour  and  fill  study.  The  model  predicts  satisfactorily 
the  evolution  of  the  water  surface  and  bed  elevations.  In  another 

test  the  model  was  used  to  simulate  bed  scour  around  a  spur  dike.  The 
shape  of  the  predicted  scour  hole  is  in  qualitative  agreement  with 
observations  reported  in  the  literature.  However,  the  deviations 
observed  in  the  predicted  locations  of  maximum  scour  and  deposition 
point  out  a  limitation  in  the  model  when  simulating  situations 

dominated  by  regions  of  flow  separation.  Further  research  is  needed 
to  correct  this  deficiency. 

5.2  RECOMMENDATIONS 

1.  It  is  recommended  that  the  model  be  further  developed  and  tested  on 
real  systems  to  ensure  its  accuracy  and  credibility.  Work  should 
continue  on  the  computer  code  to  improve  its  flexibility  and 
efficiency.  Carefully  designed  laboratory  experiments  should  be 
conducted  to  investigate  the  accuracy  and  range  of  applicability  of 
transport  algorithms  used  in  the  model  (i.e.,  turbulence  closure 
schemes,  two-dimensional  sediment  transport  functions,  etc.).  Two- 
dimensional  models  have  a  great  potential  to  study  in  detail  sediment 
related  problems  in  irregular  stream  reaches  with  significant  flow 
components  in  more  than  one  direction. 
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Verification  and  validation  of  two-dimensional  models  require  data 
with  a  degree  of  detail  and  spatial  resolution  that  is  practically 
nonexistent.  It  is  thus  recommended  that  laboratory  data  he  collected 
in  scaled  down  physical  models  reproducing  conditions  observed  in  the 
field.  Laboratory  studies  can  provide,  at  a  reasonable  cost, 
velocity,  sediment  transport,  and  cross  sectional  data  with  the  high 
degree  of  resolution  needed  in  mode]  validation.  Then,  a  few 
carefully  selected  prototype  measurements  in  the  field  will  suffice 
for  model  verification. 

It  is  recommended  that  hypothetical  situations  be  used  to  confirm  that 
the  model  responds  in  a  realistic  manner.  To  this  effect,  the  two- 
dimensional  model  can  tie  linked  to  one-dimensional  channel  and 
sediment  yield  models  to  investigate  the  dynamic  response  of  local 
bank  stabilization  structures  to  changes  in  upstream  land  management 
practices,  and  to  series  of  intense  storm  events.  These  consolidated 
models  could  be  used,  for  instance,  to  look  at  (a)  bed  and  bank 
response  in  the  vicinity  of  toe  armor,  hard  points,  fences,  etc.;  (b) 
degree  of  bank  destabilization  caused  by  proximity  of  point  bar;  etc. 
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ADDENDUM  1:  DERIVATION  OF  TWO-DIMENSIONAL  DEPTH-AVERAGED  FLOW  EQUATIONS 

The  fluid  flow  in  a  natural  stream  is  always  turbulent,  three- 
dimensional  and  time-dependent.  In  many  instances,  however,  the  flow  is 
predominantly  horizontal.  That  is,  the  flow  velocity  in  both  the 
longitudinal  and  transverse  directions  are  comparable  and  much  larger  than 
the  vertical  velocity  component.  In  such  cases,  an  expedient  approximation 
to  the  three-dimensional  flow  problem  is  to  assume  bidimensionality  along 
the  two  perpendicular  horizontal  directions,  averaging  flow  properties  in 
the  vertical  direction.  The  pressure  distribution  along  the  vertical  is 
assumed  hydrostatic,  only  time-averaged  turbulent  flow  notions  are 
considered,  and  the  effects  of  small-scale  velocity  fluctuations  are 
aggregated  into  the  shear  stress  terms. 

Continuity  equation  of  water  and  Navier-Stokes  Equation  are  written 


3u  3u2  ,  3vu  ^  3wu  ,  1  3p  1  ,^Txx  .  ^Xxy  .  ^Txz^  _  „ 
3t  +  3x  +  3y  +  3T  p  3x  ‘  p  (~3^~  3y  ~3TJ  "  0 


(1.1) 


3v  3uv  3v2 
3t  3x  3y 


+  §wv  1  §£  _  1  ,^jrx 

3z  p  3y  p  x  3x 


3t  3x 
_ZX  +  —XZ) 
3y  3z  1 


=  0 


(1.2) 


§£ 

3z 


+ 


Pg 


0 


(1.3) 


3u 

3x 


+ 


3v 

3y 


(1.4) 


where  x,  y  are  longitudinal  and  transverse  horizontal  coordinates,  z  is  the 
vertical  coordinate  measured  upwards  from  arbitrary  datum,  t  is  time,  u,  v, 
and  w  are  time-averaged  velocity  components  in  x,  y,  z  directions,  p  is 
density,  p  is  time-averaged  pressure,  g  is  acceleration  of  gravity,  and  i.. 
is  the  time-averaged  shear  stress  in  the  jth-direction  on  a  plane  normal  to 
the  ith-axies. 

The  assumed  boundary  conditions  are 

1.  the  pressure  on  the  free  surface  is  taken  as  constant. 

2.  the  free  surface  shear  stress  is  neglected. 

3.  the  fully  developed  flow  is  imposed  at  the  upstream. 

4.  no  slip,  no  seepage  on  the  channel  wall  and  bed. 

Then  they  can  be  translated  into  the  mathematical  formulations, 
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z  =  p(x,y) , 

C  1 

tl 

<  1 

It 

w  =  0 

z  =  n(x,y). 

p  =  PG 

1 

£1, 

1  H 

1 

x  |a  + 

T  = 

xx  9x 

xy  9y 

xy 

-  x  p- 

x  |n  + 

T  = 

xy  9x 

yy  3y 

yz 

TS  are  the 

surface 

wind 

on  the  bottom 
on  the  free  surface 


-s 

T 

x 

-s 

T 

y 


(1.5) 

(1.6) 

(1.7) 


Similar  expression  holds  for 


bed  shear  stresses. 

Eqn.  1 . 3  yields : 


1  a  M 

J  3^  dz  =  *  J  Pgdz 


from  which: 


p  =  pQ  +  pg(n-z) 


(1.8) 


Applying  vertical  integration,  using  the  Leibnitz's  rule,  and  Eqn. 
1.8,  the  Navier-Stoke ' s  equations  can  be  transformed  into  Eqns .  1.9-1.12 
shown  below. 

The  x-momentum  equation  becomes 

9<u>  ,  „  .  n  .  3n  ,  ,  dt  9<u2>  .  ,  .  9n 

-  u(x,y,n,t).)  gt‘  +  u  ( x ,  y ,  ,  t )  ^  +  -g— -  -  [u(x  ,y  ,r|,  t )  J 2  ^ 

,  r  r  .>|2  .  9<uv>  ,  . ,  ,  ,  9n 

+  [u(x,y,n,t))2  +  -g- -  u(x.y,n,L)  v  (x,y,n,L)  ^ 

+  u(x,y,C,t)  v  (x,y,t,t)  ||  +  u(x,y,n,t)  w  (x,y,n,t)  - 

9<t  > 

u  ( x ,  y ,  f, ,  t )  w  ( x ,  y ,  f. ,  t )  —  g~-  - 


( x ,  y ,  r) ,  t )  +  r  ( x ,  v ,  f. ,  1 1  +  -  — - 

J  9x  xx  -  9x  9y 


9<t  > 


XX 


+  T  (x.y.r.t)  > !  =  0 


t  (x,y,n,t)  +  r  (x,y,f,t) 

xy  9y  xy  9y 


9x 

> 

~9z 


(1.9) 
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3(hU)  .  3(hU2)  ,  3(hUV)  ,  w3n  .  ,  ,  3n  ,  .  3n, 

___  +  ____  +  .  u(n)[^J  +  u(n)  d  +  V(n)  dl  + 


3x 


3y 


3t 


3x 


3yJ 


u ^p)  (ft  +  +  v(p)  +  u(n)  w(n)  -  u(p)  w(P) 


3X  dy‘ 


* +  8-}^- -l  -  l  1  V>  -  V  <n)  g 

-  'xx  (r>>  IJl  *  |  |Txz<p)  -  'xy(P)  If  -  'xx  <p)  If  I  =  C 
Imposing  the  boundary  conditions,  Eqns .  1 .5-1  -  7  ,  we  have 

3(hT  ) 


3(hU)  .  3(hU2)  ,  3(hUV)  .  ,  3q  I  ouxx 

— ^ ~ -  +  — -  +  gh  r-1  -  -  — r - 


3t 


3x 


3y 


3x  p  3x 


-  ^  (Is  -  tb)  =  0 
p  X  X 


In  a  similar  fashion  the  Y-momentum  equation  yields 

3(hV)  +  3(hUV)  +  3(hV2)  +  3q  _  1  9(hTxx) 


3t 


3x 


3y 


3y  p  3x 


-  i  (ts  -  xb)  =  0 

p  y  y 


Finally,  the  continuity  equation  yields 


,  3 (hT  ) 
I  xy 

P  3y 


.  3(hT  ) 

1  _ JQL 

p  9y 


3<u>  3n  .  3C  ,  3<v>  ,  x  3n  .  ,«-x  3t 

~dT  '  u(n)  +  u( &  d  +  -JZT  -  v(n)  d  +  v<0  d 


3x 


3x  3y 


3y 


3y 


+  w(n)  -  w(£)  =  o 

3(hU)  _  3(hV) 


3x  +  "tr" +  lw(n) " u(n)  511 " 


3x  V(n)  l?J  ‘  ,W(^ 


S  -  V(D  |]  =  ° 


3(hU) 

.  3(hV) 

+  §a  . 

3x 

3y 

3t 

3(hU) 

.  3(hV) 

-  ^  = 

3x 

3y 

3t 

(1.10) 


(1.11) 


(1.12) 
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ADDENDUM  2 


COMPUTER  PROGRAM  OF  DEPTH-AVERAGED  FLOW  MODEL 
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